/*
 * Copyright 2014 Christian Bender.  All rights reserved.
 */

// standard utilities and systems includes
#include "oclDistroFix.h"


// cutoff_value
#define PROX_VALUE 0.00781250 // 1/(2^4)
#define REJECT_VALUE 512 // 128*2^4

int main(int argc, char** argv) 
{
    std::string inputfile; 
    int numthreads=8;
    inputfile="nofile";
    int mode=0; // 0=normal, 1=newcpu, 2=oldcpu
    bool benchmark = false;
    // read command line input
    // go through every command line argument
    // and see if it matches
    for(int i=0; i<argc; i++) {
        if (strcmp(argv[i], "-i")==0) {
            if ((i+1)<argc) {
                inputfile=argv[i+1];
            } else {
                std::cout << "You need to provide a Filename after the -i argument. Try oclDistroFix -i filename" << '\n';
                exit(0);  
            }
            
        }
        
        if (strcmp(argv[i], "-n")==0) {
            if ((i+1)<argc) {
                numthreads=atoi(argv[i+1]);
            } else {
                std::cout << "You need to provide the Thread number after the -n argument. Try oclDistroFix -n 8" << '\n';
                exit(0);  
            }
            
        }
        
        if (strcmp(argv[i], "-benchmark")==0) {
            benchmark=true;
        }
        
        if (strcmp(argv[i], "-oldcpu")==0) {
            mode=2;
        }
        
        if (strcmp(argv[i], "-cpu")==0) {
            mode=1;
        }
        
    }
    if (inputfile.compare("nofile")==0) {
        std::cout << "No input file name has been provided. Try oclDistroFix -i filename" << '\n';
        exit(0);
    }
    
    // open the input file
    std::ifstream file(inputfile.data());
    std::cout << "Reading Input File..." << '\n';
    // number of input lines, stored in the first line, current input line
    unsigned int inputsize, currentinputline; 
    currentinputline=0;
    
    // variables for every input line
    std::string str; 
    std::string stringarray[23];

    //read the first line
    std::getline(file, str);
    inputsize = ::atoi(str.c_str());
    
    // allocate memory
    input_data data[inputsize];
    isotope outarray[inputsize];
    
    // read the input file line by line
    while (std::getline(file, str))
    {
        if (str.at(0)=='#') {
            std::cout << "Options Input Line!" << str << '\n';
        } else {
            
            std::cout << "Reading Line " << currentinputline+1 << " of " << inputsize << '\n';
            
            int firstchar=0;
            int pos=0;
            for (unsigned int i=0; i<str.length(); i++){
                if(str[i]==',') {
                    stringarray[pos]=str.substr(firstchar,i-firstchar);
                    firstchar=i+1;
                    pos++;
                }
                stringarray[pos] = str.substr(firstchar,str.length());
                
            }
            data[currentinputline].id = ::atoi(stringarray[0].c_str());

            if (stringarray[1].compare("normal")==0) {
                data[currentinputline].mode = -1;
            } else if (stringarray[1].compare("modified")==0) {
                data[currentinputline].mode = -2;
            } else {
                data[currentinputline].mode = ::atoi(stringarray[1].c_str());
            }    
            
            data[currentinputline].C12 = ::atof(stringarray[2].c_str());
            data[currentinputline].C13= 1-data[currentinputline].C12;
            data[currentinputline].H1 = ::atof(stringarray[4].c_str());
            data[currentinputline].H2= 1-data[currentinputline].H1;
            data[currentinputline].N14 = ::atof(stringarray[6].c_str());
            data[currentinputline].N15= 1-data[currentinputline].N14;
            data[currentinputline].O16 = ::atof(stringarray[8].c_str());
            data[currentinputline].O17 = ::atof(stringarray[9].c_str());
            data[currentinputline].O18= 1-data[currentinputline].O17-data[currentinputline].O16;
            data[currentinputline].P31 = 1;
            data[currentinputline].S32 = ::atof(stringarray[12].c_str());
            data[currentinputline].S33 = ::atof(stringarray[13].c_str());
            data[currentinputline].S34 = ::atof(stringarray[14].c_str());
            data[currentinputline].S36= 1 - data[currentinputline].S32
                                          - data[currentinputline].S33
                                          - data[currentinputline].S34;
            data[currentinputline].numC = ::atoi(stringarray[16].c_str());
            data[currentinputline].numH = ::atoi(stringarray[17].c_str());
            data[currentinputline].numN = ::atoi(stringarray[18].c_str());
            data[currentinputline].numO = ::atoi(stringarray[19].c_str());
            data[currentinputline].numP = ::atoi(stringarray[20].c_str());
            data[currentinputline].numS = ::atoi(stringarray[21].c_str());
            data[currentinputline].accuracy = ::atof(stringarray[22].c_str());
            
            currentinputline++;
        }
        
        
    }
    
    std::cout << "Input File has been read!" << '\n';
    
    
    switch(mode) {
            case 0:
                runGPUImpl(data, inputsize, outarray, benchmark);
                break;
            case 1:
                runNewCPUImplParallel(data, inputsize, outarray, numthreads);
                break;
            case 2:
                runOldCPUImplParallel(data, inputsize, outarray, numthreads);
                break;
    }
 
    
    // write the output file
    std::string outputfile; 
    outputfile=inputfile+".out";
    std::ofstream outfile(outputfile.data());
    
    for (unsigned int i=0;i<inputsize;i++) {
        
        // normalize the data
        normalize_peaks(outarray[i].next);
        
        
        outfile<<data[i].id<<",";
        if (data[i].mode==-1) {
            outfile<<"normal";
        } else if (data[i].mode==-2) {
            outfile<<"modified";
        } else {
            outfile<<data[i].mode;
        }
        
        while (outarray[i].next) {
            isotope* smallest = outarray[i].next;
            //search for the smallest m/z ratio
            for (isotope* p=outarray[i].next;p;p=p->next) {
                if (p->mass<smallest->mass) smallest = p;
            }
            outfile<<","<<smallest->mass<<","<<smallest->p;
            outarray[i].next= remove_peak(outarray[i].next,smallest);
        }
        //for (isotope* p=outarray[i].next;p;p=p->next) {
        //    outfile<<","<<p->mass<<","<<p->p;
        //}
        
        outfile<<"\n";
    }
    
    outfile.flush();
    outfile.close();
}

void runNewCPUImplParallel(input_data data[], unsigned int size, isotope outarray[], int numthreads) {
    omp_set_num_threads(numthreads);
    #pragma omp parallel for
    for (unsigned int ifor = 0; ifor < size; ifor++){
        runNewCPUImpl(&data[ifor], &outarray[ifor]);
    }
}

void runGPUImpl(input_data data[], unsigned int size, isotope outarray[], bool benchmark) {
    std::chrono::nanoseconds result(0);
    oclFFT impl;
    for (unsigned int ifor = 0; ifor < size; ifor++){
        printf("Bearbeite Datensatz %i!\n",ifor);
        impl.setInput(&data[ifor]);
        //impl.calculate(&outarray[ifor]);
        result.operator +=(impl.calculate(&outarray[ifor]));
    }
    if (benchmark) {
        std::cout << "Runtime in nanoseconds: " << std::chrono::duration_cast<std::chrono::nanoseconds>(result).count() << "\n";
    }
}

void runNewCPUImpl(input_data* data, isotope* out) {
    cpuFFT impl;
    impl.setInput(data);
    out->next = impl.calculate();
}

void free_list(isotope *target)
{
    isotope* cursor;
    while(target)
    {
      cursor=target;
      target=target->next;
      free(cursor);
    }
}

void free_list(compound *target)
{
    compound* cursor;
    while(target)
    {
      cursor=target;
      target=target->next;
      free_list(cursor->isotopes);
      free(cursor);
    }
}

isotope *remove_peak(isotope *base,isotope *peak)
{
    if (!base)  {
        return NULL;
    }
    if (peak==base) {
        isotope* dummy=base;
        if(base->next)
        base->next->previous=NULL;
        base=base->next;
        free(dummy);
        return base;
    }
    peak->previous->next=peak->next;
    if (peak->next) peak->next->previous=peak->previous;
    free(peak);
    return base;
}

isotope *add_peak_sorted(isotope *base, isotope *peak) {
    // if there is no peak in base, simply return the new peak as the isotope list
    if (!base) {
        peak->previous = NULL;
        peak->next = NULL;
        return peak;
    }
    isotope* dummy = base;
    // search for a peak in the isotope list that is in close proximity to the new one
    // if a peak is found, its propability is added to the peak and it is removed from the list
    while (dummy) {
        if ((fabs(peak->mass - dummy->mass) < PROX_VALUE)) {
            peak->p += dummy->p;
            base = remove_peak(base, dummy);
            break;
        }
        dummy = dummy->next;
    }
    // the peak is now added to the list in the right position,
    // to sort the List by descending propability values
    dummy = base;
    while (dummy) {
        if (!dummy->next){
            // found the last element in the list,
            // peak will be added at its end
            dummy->next=peak;
            peak->previous=dummy;
            peak->next=NULL;
            break;
        }
            
        if (dummy->p<peak->p) {
            // found the right position,
            // insert peak in front of the current position 
            if (dummy->previous) {
                // Peak is not added as the first Element
                dummy->previous->next=peak;
            } else {
                //peak is the first element in the list
                base=peak;
            }
            peak->previous=dummy->previous;
            dummy->previous=peak;
            peak->next=dummy;
            break;
        }
        dummy = dummy->next;
    }
    return base;
}

isotope *reject_peaks(isotope* peaks)
{
    int i=0;
    
    for(isotope* p=peaks;p;p=p->next){
        if(++i>=REJECT_VALUE) {
                free_list(p->next);
                p->next=NULL;
        }
    }
    return peaks;
}

void runOldCPUImplParallel(input_data data[], unsigned int size, isotope outarray[], int numthreads) {
    omp_set_num_threads(numthreads);
    #pragma omp parallel for
    for (unsigned int ifor = 0; ifor < size; ifor++){
        runOldCPUImpl(&data[ifor], &outarray[ifor]);
    }
}

void count_peaks(isotope* peaks) {
    int i=0;
    for (isotope* p=peaks;p;p=p->next) {
        i++;
    }
    printf("Es wurden %i Datensätze gefunden!\n",i);
}

void normalize_peaks(isotope* peaks) {
    if (peaks) {
    cl_double factor = (100.0/peaks->p);
        for (isotope* p=peaks;p;p=p->next) {
            p->p=p->p*factor;
    }
    }
}

void runOldCPUImpl(input_data* data, isotope* out) {
    
        isotope *ci; // current isotope being added
        isotope *icur; // cursor to the last element in the list
        compound *cc; //current compound being added
        compound *ccur; // cursor to the last compound in the list
        
        // compound creation
        compound *verbindung=NULL;
        cc = (compound*) malloc(sizeof(compound));
        cc->amount=data->numC;
        cc->previous=NULL;
        cc->next=NULL;
        verbindung = cc;

        // c12
        ci = (isotope*)malloc(sizeof(isotope));
        ci->p = data->C12;
        ci->mass = C12_MASS;
        ci->previous=NULL;
        ci->next=NULL;
        cc->isotopes=ci;
        icur=ci;
        
        // c13
        ci = (isotope*)malloc(sizeof(isotope));
        ci->p = data->C13;
        ci->mass = C13_MASS;
        ci->previous=icur;
        icur->next=ci;
        ci->next=NULL;
        icur=ci;
        
        // helium
        ccur = cc;
        cc=(compound*)malloc(sizeof(compound));
        ccur->next=cc;
        cc->previous=ccur;
        cc->amount = data->numH;
        cc->next=NULL;
        
        // h1
        ci = (isotope*)malloc(sizeof(isotope));
        ci->p = data->H1;
        ci->mass = H1_MASS;
        ci->previous=NULL;
        cc->isotopes=ci;
        icur=ci;
        
        // h2
        ci = (isotope*)malloc(sizeof(isotope));
        ci->p = data->H2;
        ci->mass = H2_MASS;
        ci->previous=icur;
        icur->next=ci;
        ci->next=NULL;
        icur=ci;
        
        // nitrogen
        ccur = cc;
        cc=(compound*)malloc(sizeof(compound));
        ccur->next=cc;
        cc->previous=ccur;
        cc->amount=data->numN;
        cc->next=NULL;
        
        // n14
        ci = (isotope*)malloc(sizeof(isotope));
        ci->p = data->N14;
        ci->mass = N14_MASS;
        ci->previous=NULL;
        ci->next=NULL;
        cc->isotopes=ci;
        icur=ci;
        
        // n15
        ci = (isotope*)malloc(sizeof(isotope));
        ci->p = data->N15;
        ci->mass = N15_MASS;
        ci->previous=icur;
        icur->next=ci;
        ci->next=NULL;
        icur=ci;
        
        // oxygen
        ccur = cc;
        cc=(compound*)malloc(sizeof(compound));
        ccur->next=cc;
        cc->previous=ccur;
        cc->amount = data->numO;
        cc->next=NULL;
        
        // 016
        ci = (isotope*)malloc(sizeof(isotope));
        ci->p = data->O16;
        ci->mass = O16_MASS;
        ci->previous=NULL;
        cc->isotopes=ci;
        icur=ci;
        
        // o17
        ci = (isotope*)malloc(sizeof(isotope));
        ci->p = data->O17;
        ci->mass = O17_MASS;
        ci->previous=icur;
        icur->next=ci;
        ci->next=NULL;
        icur=ci;
        
        // o18
        ci = (isotope*)malloc(sizeof(isotope));
        ci->p = data->O18;
        ci->mass = O18_MASS;
        ci->previous=icur;
        icur->next=ci;
        ci->next=NULL;
        icur=ci;
        
        // phosporus
        ccur = cc;
        cc=(compound*)malloc(sizeof(compound));
        ccur->next=cc;
        cc->previous=ccur;
        cc->amount = data->numP;
        cc->next=NULL;
        
        // p31
        ci = (isotope*)malloc(sizeof(isotope));
        ci->p = data->P31;
        ci->mass = P31_MASS;
        ci->previous=NULL;
        ci->next=NULL;
        cc->isotopes=ci;
        icur=ci;
        
        // sulfur
        ccur = cc;
        cc=(compound*)malloc(sizeof(compound));
        ccur->next=cc;
        cc->previous=ccur;
        cc->amount = data->numS;
        cc->next=NULL;
        
        // s32
        ci = (isotope*)malloc(sizeof(isotope));
        ci->p = data->S32;
        ci->mass = S32_MASS;
        ci->previous=NULL;
        ci->next=NULL;
        cc->isotopes=ci;
        icur=ci;
        
        // s33
        ci = (isotope*)malloc(sizeof(isotope));
        ci->p = data->S33;
        ci->mass = S33_MASS;
        ci->previous=icur;
        icur->next=ci;
        ci->next=NULL;
        icur=ci;
        
        // s34
        ci = (isotope*)malloc(sizeof(isotope));
        ci->p = data->S34;
        ci->mass = S34_MASS;
        ci->previous=icur;
        icur->next=ci;
        ci->next=NULL;
        icur=ci;
        
        // s36
        ci = (isotope*)malloc(sizeof(isotope));
        ci->p = data->S36;
        ci->mass = S36_MASS;
        ci->previous=icur;
        icur->next=ci;
        ci->next=NULL;
        icur=ci;
        
        // the resulting peaks
        isotope *peaks;

        // Calculate Peaks
        compound *c;
        isotope *npeaks, *p, *i, *np1;

        // initializing List
        peaks = (isotope*) malloc(sizeof (isotope));
        peaks->mass = 0;
        peaks->p = 1;
        peaks->previous = NULL;
        peaks->next = NULL;
        
        printf("Data Structure created successfully. ID: %i  \n",data->id);

        for (c = verbindung; c; c = c->next) {
            for (int anzahl = 0; anzahl < c->amount; anzahl++) {
                npeaks = NULL;

                for (p = peaks; p; p = p->next) {
                    for (i = c->isotopes; i; i = i->next) {
                        np1 = (isotope*) malloc(sizeof (isotope));
                        np1->mass = p->mass + i->mass;
                        np1->p = p->p * i->p;
                        np1->next=NULL;
                        np1->previous=NULL;
                        npeaks = add_peak_sorted(npeaks, np1);
                        np1=NULL;
                        
                    }
                }
                free_list(peaks);
                peaks = npeaks;
                peaks = reject_peaks(peaks);
            }
        }
        
        out->next=peaks;
        free_list(verbindung);
        peaks=NULL;

}
